source('../init.R')

## plot the density function of wilcoxon U
## shade in 2-sided tails given by alpha
my.plot.wilcox <- function(n1, n2, alpha, ...) {
    .max.u <- qwilcox(1, n1, n2)
    .seq.u <- 0:.max.u
    ## max density, at p=50%
    .max.dens <-  dwilcox(qwilcox(0.5, n1, n2), n1, n2) 
    ## tail position
    .lower <- qwilcox(alpha/2, n1, n2)
    .upper <- qwilcox(1-alpha/2, n1, n2)
    .df = data.frame(U=.seq.u, Density=dwilcox(.seq.u, n1, n2))
    xyplot(Density ~U, .df, type=c('l', 'g'), col='black', 
        panel=function(x,y,...) { 
            ## vertical lines at alpha/2
            panel.abline(v=.lower, col='red'); 
            panel.abline(v=.upper, col='red'); 
            panel.xyplot(x,y,...)
            ## shade tails
            panel.xyplot(x[x < .lower],y[x < .lower],type='h') 
            panel.xyplot(x[x > .upper],y[x > .upper],type='h') 
        }, 
        ylim=c(0, .max.dens),
        ...
    )
}



## for each n and alpha, make a plot
.n=c(10, 20)
.alpha=c(0.05, 0.01)
l_ply(.alpha, function(..a) {
    l_ply(.n, function(..n) {
        .file <- sprintf('fig-wk7-wilcox-n%d-alpha%0.2f.png', ..n, ..a)
        mypng(.file,  code=
            my.plot.wilcox(..n, ..n, ..a, 
                sub=sprintf('Sampling distribution of Wilcoxon U for N=%d \n Showing 2-sided tails for alpha=%0.2f', ..n, ..a)
            )
        )
    })
})
